######
#This portion of the code simply merges ideal point and alliance data onto our events data#
######

rm(list=ls())
setwd("~/Dropbox/Foreign Gifts/Replication materials")
library(devtools)
#The version of record of this paper uses the countrycode package as of the git commit dated April 2, 2015 and found at https://github.com/vincentarelbundock/countrycode/commit/893d2319f542e411ab70f117ae5d25709747012e  Earlier versions of the countrycode package contain errors that are consequential for merging the data.  It is conceivable that later versions will also impact merging behavior, although we believe that this version is accurate for the purposes of the operations we conduct.  Nonetheless, replicators should install that version from github to reproduce the results here.
install_github('vincentarelbundock/countrycode',ref="893d2319f542e411ab70f117ae5d25709747012e")
library(countrycode)

load('internationalEvents.RData')


#We use Version 7.1 of the Bailey, Strezhnev, and Voeten dyadic ideal point data (date May 11, 2015).  Replicators may independently obtain this file at https://dataverse.harvard.edu/dataset.xhtml?persistentId=hdl:1902.1/12379&version=7.1  This version of the data must be used  to reproduce the results here.  The file loaded here is an unmodified version of the original.

load('Idealpointsdyadicdistance_1.RData')

nonUSData<-internationalEvents[internationalEvents[,"actor1country"]!="UNITED STATES"&internationalEvents[,"actor2country"]!="UNITED STATES",]
nonUSData[,"actor1country"]<-countrycode(nonUSData[,"actor1country"],"country.name","country.name")
nonUSData[,"actor2country"]<-countrycode(nonUSData[,"actor2country"],"country.name","country.name")

idealPoints<-x[x[,"countryname1"]=="United States of America",c("year","cowcode2","idealpointdistance")]
idealPoints<-data.frame(idealPoints,"my"=floor((idealPoints[,"year"]-1)/4))
aggDist<-aggregate(idealpointdistance~my+cowcode2,data=idealPoints,FUN=mean)
colnames(aggDist)[3]<-"aggDist"
idealPoints<-merge(idealPoints,aggDist)
idealPoints<-idealPoints[,c("year","cowcode2","aggDist")]
colnames(idealPoints)<-c("year","actor1country","idealpointdistance.1")
idealPoints[,"actor1country"]<-countrycode(idealPoints[,"actor1country"],"cown","country.name")
nonUSData<-merge(nonUSData,idealPoints,all.x=T)
colnames(idealPoints)<-c("year","actor2country","idealpointdistance.2")
nonUSData<-merge(nonUSData,idealPoints,all.x=T)

#This is an unmodified copy of version 4.1 of the COW alliances dataset.
alliances<-read.csv("COWAlliances.csv")
#Note that the US appears first in all alliances of which it is a member in this data
usalliances<-alliances[(alliances[,"ccode1"]==2&alliances[,"year"]>1945),]
usalliances[,"state_name2"]<-countrycode(usalliances[,"ccode2"],"cown","country.name")
usalliances<-usalliances[,c("state_name2","year","defense")]
usalliances<-aggregate(defense~state_name2+year,data=usalliances,FUN=sum)
colnames(usalliances)<-c("actor1country","year","ally.1")

#Note to replicators: This merge results in a "1" for ally.1 and ally.2 whenever a country IS allied with the United States, but often produces an "NA" when a country is NOT allied with the United States.  An alternative merge should be used to run an analysis on countries NOT allied with the US (no such results are presented in the paper or supplementary information)
nonUSData<-merge(nonUSData,usalliances,all.x=T)
colnames(usalliances)<-c("actor2country","year","ally.2")
nonUSData<-merge(nonUSData,usalliances,all.x=T)

USData<-internationalEvents[internationalEvents[,"actor1country"]=="UNITED STATES"|internationalEvents[,"actor2country"]=="UNITED STATES",]


#####
##Begining here we replicate the results of the main paper##
#####


library(parallel)
library(lubridate)
library(xts)
library(nlme)
library(lmtest)
library(tseries)
library(sandwich)


cutoff<-median(c(nonUSData[complete.cases(nonUSData[,c("idealpointdistance.1","idealpointdistance.2")]),"idealpointdistance.1"],nonUSData[complete.cases(nonUSData[,c("idealpointdistance.1","idealpointdistance.2")]),"idealpointdistance.2"]))
inUSBloc<-nonUSData[nonUSData[,"idealpointdistance.1"]<=cutoff&nonUSData[,"idealpointdistance.2"]<=cutoff&!is.na(nonUSData[,"idealpointdistance.1"])&!is.na(nonUSData[,"idealpointdistance.2"]),]
##
usBlocTS<-xts(inUSBloc[,"GOLDSTEINSCALE"],
              order.by=as.Date(paste(inUSBloc[,"year"], substr(inUSBloc$date, 6, 7), "01", sep = "-"))
                  )
monthlyMean<-apply.monthly(usBlocTS,mean)
monthlyCount<-apply.monthly(usBlocTS,nrow)
isElection<-year(monthlyCount)%%4==0
inWindow<-month(monthlyCount)%in%c(6,7,8,9,10,11)
monthlyCount<-data.frame(monthlyCount,isElection)

#Total count as reported in text of ally-ally interactions

#Year level count as reported in text 
t.test(monthlyCount~isElection,data=monthlyCount)
#"Within window" count as reported in text
t.test(monthlyCount~isElection,data=monthlyCount[inWindow,])

#T-test at annual level as reported in text 
t.test(monthlyMean~isElection)
#T-test in window as reported in text
t.test(monthlyMean[inWindow]~isElection[inWindow])

#Augmented Dickey Fuller as reported in text
adf.test(monthlyMean)

proximateElection<-isElection&inWindow
proximateMidterm<-year(monthlyMean)%%4==2&inWindow
theMonth<-month(monthlyMean)
theYear<-year(monthlyMean)
augmented<-cbind("meanScale"=monthlyMean,"proximate"=proximateElection,"premidterm"=proximateMidterm,"Month"=theMonth,"Year"=theYear)

augmented<-cbind(augmented,NA)
#Fixed effects are assembled at the President-term level - i.e., switch to a new FE
#When President changes (incl. via death) or after a subsequent innauguration
colnames(augmented)[ncol(augmented)]<-"Term"
augmented["1933-03-01/1937-01-31","Term"]<-1
augmented["1937-02-01/1941-01-31","Term"]<-2
augmented["1941-02-01/1945-03-31","Term"]<-3
augmented["1945-04-01/1949-01-31","Term"]<-4
augmented["1949-02-01/1953-01-31","Term"]<-5
augmented["1953-02-01/1957-01-31","Term"]<-6
augmented["1957-02-01/1961-01-31","Term"]<-7
augmented["1961-02-01/1963-11-31","Term"]<-8
augmented["1963-12-01/1965-01-31","Term"]<-9
augmented["1965-02-01/1969-01-31","Term"]<-10
augmented["1969-02-01/1973-01-31","Term"]<-11
augmented["1973-02-01/1974-08-01","Term"]<-12
augmented["1974-08-01/1977-01-31","Term"]<-13
augmented["1977-02-01/1981-01-31","Term"]<-14
augmented["1981-02-01/1985-01-31","Term"]<-15
augmented["1985-02-01/1989-01-31","Term"]<-16
augmented["1989-02-01/1993-01-31","Term"]<-17


#########
#Replication of 6 Models in Table 3
########

model1<-lm(meanScale~proximate+as.character(Month),data=augmented)
coeftest(model1,vcov = vcovHAC(model1))
#
model2<-lm(meanScale~proximate+as.character(Month)+as.character(ceiling(Year/2)),data=augmented)
coeftest(model2,vcov = vcovHAC(model2))
#
model3<-lm(meanScale~proximate+as.character(Month)+as.character(floor((Year-1)/4)),data=augmented)
coeftest(model3,vcov = vcovHAC(model3))

#
model4<-lm(meanScale~proximate+premidterm+as.character(Month),data=augmented)
coeftest(model4,vcov = vcovHAC(model4))
#
model5<-lm(meanScale~proximate+premidterm+as.character(Month)+as.character(ceiling(Year/2)),data=augmented)
coeftest(model5,vcov = vcovHAC(model5))
#
model6<-lm(meanScale~proximate+premidterm+as.character(Month)+as.character(floor((Year-1)/4)),data=augmented)
coeftest(model6,vcov = vcovHAC(model6))


#Percentile comparison discussed in text
monthlyEc<-ecdf(as.numeric(augmented[,"meanScale"]))
monthlyEc(median(augmented[,"meanScale"])-0.25)


crossBloc<-nonUSData[(!(nonUSData[,"idealpointdistance.1"]<=cutoff&nonUSData[,"idealpointdistance.2"]<=cutoff)&(nonUSData[,"idealpointdistance.1"]<=cutoff|nonUSData[,"idealpointdistance.2"]<=cutoff))&(!is.na(nonUSData[,"idealpointdistance.1"])&!is.na(nonUSData[,"idealpointdistance.2"])),]
#Number of cross bloc interactions as reported in text
nrow(crossBloc)

crossTS<-xts(crossBloc[,"GOLDSTEINSCALE"],order.by=crossBloc[,"date"])
crossMean<-apply.monthly(crossTS,mean)
crossCount<-apply.monthly(crossTS,nrow)
isElection<-year(crossCount)%%4==0
inWindow<-month(crossCount)%in%c(6,7,8,9,10,11)
crossCount<-data.frame(crossCount,isElection)

#Comparison of event count (not reported in text)
t.test(crossCount~isElection,data=crossCount)
t.test(crossCount~isElection,data=crossCount[inWindow,])

proximateElection<-isElection&inWindow
proximateMidterm<-year(crossMean)%%4==2&inWindow
theMonth<-month(crossMean)
theYear<-year(crossMean)
augmentedCross<-cbind("meanScale"=crossMean,"proximate"=proximateElection,"premidterm"=proximateMidterm,"Month"=theMonth,"Year"=theYear)

## Table 4
crossMod1<-lm(meanScale~proximate+as.character(Month),data=augmentedCross)
coeftest(crossMod1,vcov = vcovHAC(crossMod1))
crossMod2<-lm(meanScale~proximate+as.character(Month)+as.character(floor((Year-1)/4)),data=augmentedCross)
coeftest(crossMod2,vcov = vcovHAC(crossMod2))
crossMod3<-lm(meanScale~proximate+as.character(Month)+as.character(ceiling(Year/2)),data=augmentedCross)
coeftest(crossMod3,vcov = vcovHAC(crossMod3))

usTS<-xts(USData[,"GOLDSTEINSCALE"],order.by=USData[,"date"])
usMean<-apply.monthly(usTS,mean)
usCount<-apply.monthly(usTS,nrow)
isElection<-year(usCount)%%4==0
inWindow<-month(usCount)%in%c(6,7,8,9,10,11)
theYear<-year(usCount)
theMonth<-month(usCount)
proximateMidterm<-inWindow&theYear%%4==2
usCount<-data.frame(usCount,isElection)
#Unreported count comparison - NOTE: number of US events does vary in election years
t.test(usCount~isElection,data=usCount)
t.test(usCount~isElection,data=usCount[inWindow,])
augmentedUS<-cbind("meanScale"=usMean,"proximate"=isElection&inWindow,"premidterm"=proximateMidterm,"Month"=theMonth,"Year"=theYear)

## Table 4
usMod1<-lm(meanScale~proximate+as.character(Month),data=augmentedUS)
coeftest(usMod1,vcov = vcovHAC(usMod1))
usMod2<-lm(meanScale~proximate+as.character(Month)+as.character(ceiling(Year/2)),data=augmentedUS)
coeftest(usMod2,vcov = vcovHAC(usMod2))
usMod3<-lm(meanScale~proximate+as.character(Month)+as.character(floor((Year-1)/4)),data=augmentedUS)
coeftest(usMod3,vcov = vcovHAC(usMod3))
